/*
 * ExactDiag.cpp
 *
 *  Created on: Oct 9, 2012
 *      Author: snirgaz
 */

#include "ExactDiag.h"

vector<int> mzStates(int L, int mz) {
	vector<int> mzStates;
	boost::dynamic_bitset<> s(L);
	for (int n = 0; n < L; n++)
		s[n] = 1;
	unsigned int maxState = s.to_ulong();
	for (unsigned int n = 0; n <= maxState; n++) {
		boost::dynamic_bitset<> s(L, n);
		int mzc = 0;
		for (int site = 0; site < L; site++)
			mzc += s[site];
		if (mzc == mz)
			mzStates.push_back(n);
	}
	return mzStates;
}

void ExactDiag::init(int L, double alpha) {
	L_ = L;
	alpha_ = alpha;
	outPath_ = "/home/snirgaz/sim_data/exact_diag/"
			+ boost::lexical_cast<string>(L) + "_"
			+ boost::lexical_cast<string>(alpha_);
	boost::filesystem::create_directories(outPath_);
}

void ExactDiag::constructH() {
	numOfmzStates_ = mzStates_.size();
	H_ = MatrixXd::Zero(numOfmzStates_, numOfmzStates_);
	vector<int>::iterator it_mz = mzStates_.begin();
	int l = (int) sqrt(L_);
	for (; it_mz != mzStates_.end(); it_mz++) {
		int state1 = it_mz - mzStates_.begin();
		boost::dynamic_bitset<> s(L_, *it_mz);
		for (int n = 0; n < L_; n++) {
			int ixn, iyn, ixnn, iynn;
			ixn = n % l;
			iyn = n / l;
			array<int, 2> nns;
			// iy Fixed
			iynn = iyn;
			ixnn = ((ixn + 1) == l) ? 0 : (ixn + 1);
			nns[0] = ixnn + l * iynn;
			// ix Fixed
			ixnn = ixn;
			iynn = ((iyn + 1) == l) ? 0 : (iyn + 1);
			nns[1] = ixnn + l * iynn;
			for (int i = 0; i < 2; i++) {
				int nn = nns[i];
				double a = ((i == 0) && (ixn % 2)) ? alpha_ : 1.;
				if (s[n] == s[nn]) {
					H_(state1, state1) += 0.25 * a;
				} else {
					H_(state1, state1) -= 0.25 * a;
					boost::dynamic_bitset<> mask(L_);
					mask[n] = 1;
					mask[nn] = 1;
					boost::dynamic_bitset<> sflip = s ^ mask;
					unsigned int sflipint = sflip.to_ulong();
					int state2 = std::find(mzStates_.begin(), mzStates_.end(),
							sflipint) - mzStates_.begin();
					H_(state1, state2) = 0.5 * a;
				}
			}
		}
	}
}

void ExactDiag::diag() {
	SelfAdjointEigenSolver<MatrixXd> eigensolver(H_);
	string fileName = outPath_ + "/eigsys" + boost::lexical_cast<string>(mz_)
			+ ".h5";
	hdf5RW h5File;
	h5File.initNew(fileName);
	DimsHdf5 matSize(2, eigensolver.eigenvectors().rows());
	h5File.writeComp("/vec", matSize, eigensolver.eigenvectors().data());
	DimsHdf5 vecSize(1, eigensolver.eigenvalues().size());
	h5File.writeComp("/val", vecSize, eigensolver.eigenvalues().data());
	cout << "mz =" << mz_ << " num Of States " << numOfmzStates_ << " E - "
			<< eigensolver.eigenvalues().minCoeff() << "\n";

}
void ExactDiag::diagAll() {
	for (mz_ = 0; mz_ <= L_; mz_++) {
		mzStates_ = mzStates(L_, mz_);
		constructH();
		diag();
	}
}

void LFCorr::init(int L, double beta, int maxOmega, double alpha) {
	L_ = L;
	beta_ = beta;
	outPath_ = "/home/snirgaz/sim_data/exact_diag/"
			+ boost::lexical_cast<string>(L) + "_"
			+ boost::lexical_cast<string>(alpha);
	boost::filesystem::create_directories(outPath_);
	maxOmega_ = maxOmega;
	corr_ = VectorXd::Zero(maxOmega_);

}
void LFCorr::readData(int mz) {
	numOfmzStates_ = mzStates_.size();
	string fileName = outPath_ + "/eigsys" + boost::lexical_cast<string>(mz_)
			+ ".h5";
	hdf5RW h5File;
	h5File.initRead(fileName);
	eigensys_.eigVal.resize(numOfmzStates_);
	eigensys_.eigVec.resize(numOfmzStates_, numOfmzStates_);
	h5File.read("/vec", eigensys_.eigVec.data());
	h5File.read("/val", eigensys_.eigVal.data());
	cout << "mz =" << mz_ << " num Of States " << numOfmzStates_ << " E - "
			<< eigensys_.eigVal.minCoeff() << "\n";

}
void LFCorr::constructHLF() {
	HLF_ = MatrixXd::Zero(numOfmzStates_, numOfmzStates_);
	vector<int>::iterator it_mz = mzStates_.begin();
	int l = (int) sqrt(L_);
	for (; it_mz != mzStates_.end(); it_mz++) {
		int state1 = it_mz - mzStates_.begin();
		boost::dynamic_bitset<> s(L_, *it_mz);
		for (int n = 0; n < L_; n++) {
			int ixn, iyn, ixnn, iynn;
			ixn = n % l;
			iyn = n / l;
			array<int, 2> nns;
			// iy Fixed
			iynn = iyn;
			ixnn = ((ixn + 1) == l) ? 0 : (ixn + 1);
			nns[0] = ixnn + l * iynn;
			// ix Fixed
			ixnn = ixn;
			iynn = ((iyn + 1) == l) ? 0 : (iyn + 1);
			nns[1] = ixnn + l * iynn;
			for (int i = 0; i < 2; i++) {
				int nn = nns[i];
				double sign = (i == 0) ? 1 : -1;
				if (s[n] == s[nn]) {
					HLF_(state1, state1) += 0.25 * sign;
				} else {
					HLF_(state1, state1) -= 0.25 * sign;
					boost::dynamic_bitset<> mask(L_);
					mask[n] = 1;
					mask[nn] = 1;
					boost::dynamic_bitset<> sflip = s ^ mask;
					unsigned int sflipint = sflip.to_ulong();
					int state2 = std::find(mzStates_.begin(), mzStates_.end(),
							sflipint) - mzStates_.begin();
					HLF_(state1, state2) = 0.5 * sign;
					HLF_(state2, state1) = 0.5 * sign;
				}
			}
		}
	}
}
void LFCorr::CalcHLFCorr(bool sym) {
	double facSym = (sym) ? 2. : 1.;
	EigenSys &eigenSys = eigensys_;
	VectorXd eig_val_shift = eigenSys.eigVal.array() + 4.;
	int numOfmzStates = eigenSys.eigVal.size();
	MatrixXd overlap = eigenSys.eigVec.transpose() * HLF_ * eigenSys.eigVec;
	for (int m = 0; m < numOfmzStates_; m++) {
		Z_ += facSym * exp(-beta_ * eig_val_shift(m));
		for (int n = 0; n < m; n++) {
			double me = overlap(m, n) * overlap(m, n);
			double de = eig_val_shift(m) - eig_val_shift(n);
			double fac = me
					* (exp(-beta_ * eig_val_shift(m))
							- exp(-beta_ * eig_val_shift(n))) * (-2 * de);
			for (int omegam = 0; omegam < maxOmega_; omegam++) {
				double fm = (omegam * 2. * boost::math::constants::pi<double>())
						/ beta_;
				corr_(omegam) += facSym * fac / (fm * fm + de * de);
			}
		}
	}
}
void LFCorr::calcAll() {
	for (mz_ = 0; mz_ <= L_ / 2; mz_++) {
		mzStates_ = mzStates(L_, mz_);
		readData(mz_);
		constructHLF();
		bool sym = (mz_ == L_ / 2) ? false : true;
		CalcHLFCorr(sym);
	}
	string fileName = outPath_ + "/CorrHLF_2.h5";
	hdf5RW h5File;
	h5File.initNew(fileName);
	corr_ /= Z_;
	DimsHdf5 vecSize(1, corr_.size());
	h5File.writeComp("/corr", vecSize, corr_.data());
}
void EnergyMeasure::init(int L, double beta, double alpha) {
	L_ = L;
	beta_ = beta;

	outPath_ = "/home/snirgaz/sim_data/exact_diag/"
			+ boost::lexical_cast<string>(L) + "_"
			+ boost::lexical_cast<string>(alpha);
	boost::filesystem::create_directories(outPath_);

	energy_ = 0.;

}
void EnergyMeasure::readData(int mz) {
	numOfmzStates_ = mzStates_.size();
	string fileName = outPath_ + "/eigsys" + boost::lexical_cast<string>(mz_)
			+ ".h5";
	hdf5RW h5File;
	h5File.initRead(fileName);
	eigensys_.eigVal.resize(numOfmzStates_);
	eigensys_.eigVec.resize(numOfmzStates_, numOfmzStates_);
	h5File.read("/vec", eigensys_.eigVec.data());
	h5File.read("/val", eigensys_.eigVal.data());
	cout << "mz =" << mz_ << " num Of States " << numOfmzStates_ << " E - "
			<< eigensys_.eigVal.minCoeff() << "\n";

}

void EnergyMeasure::CalcEnergy(bool sym) {
	double facSym = (sym) ? 2. : 1.;
	EigenSys &eigenSys = eigensys_;
	VectorXd eig_val_shift = eigenSys.eigVal.array() + 12.;
	VectorXd boltz = (-beta_ * eig_val_shift).array().exp();
	energy_ += facSym * (boltz.dot(eig_val_shift));
	Z_ += facSym * (boltz).sum();
}

void EnergyMeasure::calcAll() {
	for (mz_ = 0; mz_ <= L_ / 2; mz_++) {
		mzStates_ = mzStates(L_, mz_);
		readData(mz_);
		bool sym = (mz_ == L_ / 2) ? false : true;
		CalcEnergy(sym);
	}
	string fileName = outPath_ + "/energy.h5";
	hdf5RW h5File;
	h5File.initNew(fileName);
	energy_ /= Z_;
	cout << energy_ - 12;
	h5File.writeSingle("/energy", energy_ - 12);
}
